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1. Introduction 



Stochastic linear programs can be formulated for a variety of applications. Some examples include airline 
scheduling (Ferguson and Dantzig [1956]), financial planning (Kusy and Ziemba [1986]), energy modeling 
(Birge [1987]) and water resource planning (Pr^kopa and Szantai [1978]). The basic model we consider here 
is the stochastic linear program with recourse in the following general form: 

min x {c r x + Q(x)\Ax = 6,x > 0} 

where 

Q(x) = I Q{x,t,4)P{dt,d*) 

and the recourse function is defined as 

Q(x, Z, <t>) = min{<7 r y|Wt/ = £ - Tx, u + <f> > y > 0}, 

where x 6 3R ni , y 6 9?"*, b 6 and (£, 4>) is a random vector on the probability space (3R m3+n3 , P) 

with support, 5x$. The vectors, c, <7, and u, and matrices, A , W , and T are dimensioned correspondingly. 
The fundamental problem in stochastic programming is to evaluate the integral of Q. In this paper, we 
describe a method for finding an upper bound on Q that requires a polynomial number of operations in the 
number of random variables. 

Previous results in bounding expressions for Q are described in Birge and Wets [1986a]. The bounds are 
based on the convexity and positive homogeneity of Q . The first result is due to Jensen [1906] *s inequality 
which provides a lower bound on Q. The usefulness of this lower bound is that it requires an evaluation of Q 
at one point (the mean of the random variables) and has been found to be generally sharp in some practical 
examples (see, e.g., Hausch and Ziemba [1983]). Madansky [1959] provided an upper bound following 
Edmundson [1956] that is based on the theory of moment spaces and amounts to weighting the extreme 
points of the support of the random variables. Ben-Tal and Hochman [1972] and Huang, Ziemba and Ben-Tal 
[1977] refined this bound for independent random variables. Dupa£ovd [1976] formulated a bound of the same 
general type for dependent random variables that was extended to unbounded ranges and non-polyhedral 
sets in Gassmann and Ziemba [1986]. FVauendorfer [1986] provided a sharper bound in the bounded range, 
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dependent variable case, and Birge and Wallace [1986] gave a bound and method for refinement for special 
cases of dependent random variables. 

The upper bounds mentioned above all have the property that they are solutions to moment problems 
with varying conditions. Dupa£ova’s work on minimax solutions (2d£kova [1966]) led to these conclusions 
and to the use of the generalized moment problem. Ermoliev, et al. [1987] provided a general programming 
framework for solving the general problem. It is used in Birge and Wets [1987] for bounds with piecewise 
linear approximations on moment constraints and in Cipra [1985] with first and second moment constraints. 

The problem with each of these bounds is that they require an exponentially increasing number of 
function evaluations as the number of random variables increases. An alternative for this situation was given 
by the ray approximation procedure in Birge and Wets [1986a]. This uses the sublinearity property of the 
recourse function to obtain a separable function that majorizes Q . This approach is generalized in Birge 
and Wets [1986b]. Wallace [1987b], on the other hand, formulated a procedure that applies to problems in 
which the recourse function involves the solution of a network problem. Our procedure is a combination and 
generalization of these two basic approaches. The algorithm we give provides a separable piecewise linear 
function that bounds Q throughout the support of the random variables and can be easily evaluated. 

Section 2 presents our basic algorithm and the separable piecewise linear upper bound ( SPLU). Its 
properties are described in Section 3. Section 4 gives an illustrative small example and provides comparison 
with the upper bound of Edmundson and Madansky. Extensions of the basic algorithm and conclusions are 
given in Section 5. 

2. The Basic Algorithm 

We give a general method for finding an upper bound on the expected value of the value of a linear 

program with random right-hand sides and random upper bounds on the variables. To simplify notation 

and to establish general results, we consider the following system : 

A\ x = &i + £ 

A 2 x = b 2 (1) 

0 < x < c + 4> 
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where A\ € 9ft miXn , A 2 € 3f^( m “ m i) Xrl J [A\\A 2 ) T = A is the coefficient matrix, (&i|&2) r = b is the fixed 
part of the right-hand side, c is the fixed part of the bounds on the variables, £ is the random availability 
of resources and <j> is the random part of the variable capacities, where <j> > 0. We assume that there is a 
positive probability that <f> = 0. Next define Q(£,<j>) by 

Q(Z,<j>) = min{^ T x|(l)}. (2) 

Finally define x(£> as the set of x-vectors satisfying 

A x x = £ 

A 2 x = 0 ( 3 ) 

dT < x < <}> + d + . 

Our goal is fo find an upper bound on Q(£,<£), or, more precisely, on EQ(£, <j>). We do this by finding a 
separable piecewise linear function {/(£, <j>) defined by 

mi ( q T x' + (ti - £)if& > & 

U{£y <f>) = Q(Z, o) + H(<j>) + < 

»=1 { q T x'-{Zi - £,)if£ < 6 

where £ t - = E£*, and H {(j>) is a piecewise linear function in <j>. 

Algorithm 1 

Step 0: Find Q(£, 0) with optimal solution x°, where 

! ej SJ* 1 (6 + £) if i is basic, 

0 if i is nonbasic at lower bound, 

c(t) if t is nonbasic at upper bound. 

Assume for simplicity that the first m variables are basic. Let x i+ = (Bq 1 e t -, 0,0, . . . ,0) and x'~ = 
(-So le »,0,0,...,0) where t = 1,2 Let 

mi mi 

a x (i) = max - z°(i) - ^ r , + (t)y J+ - 

3=2 3=2 

subject to 

yt+-y>- 

<T n <e ? - <e? iax ,y = 2,... ) m 1 . 
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Let 



subject to 


m i m i 

^ 1 (t) = min — x°(t) — ^ x 3+ (i)y ,+ — x } ~(i)y J ~ + c(i) 

3 = 2 3= 2 


for all i = 1, . . . , n. 


y J+ - y*~ = £y - 1; 

£f in < £/ < £“ ax , j = 2,...,m u 



If a 1 ^) > 0 for some i or P l [i) < 0 for some z, let x l + = ct 1 (i) = (0, . ..,0), and = z°(z) + c(z) for all 

i — 1, . . . , mi and go to Step 1 with r = 1. 

Otherwise, check 





e[i) = max — x 1+ (z)y + — x 1 (i)y 


subject to 


1 

<0 

II 

1 

+ 


and 


^r in < 6 < er x 



f(i) = min - z 1+ (z)y + - x 1 ( z)y 



subject to 


y + - y = £i - lx 




£r in < ii < er ax - 



If a 1 -f e < 0 and /9 1 + / > 0, then Q(£, <f>) is linear in £, go to Step 4. 

Otherwise, let r = 1 and go to Step 1. 

Step 1: If £™ ax < -boo, solve 

min{< 7 T a; | x((6”“ c “ | r )e r , 0, a r , &]} = g 7 V + (£ nax - | r ). 
Else (let = oo if /? r (z) = +oo, /?*(z) = 0 otherwise) and solve 



If (™ in > — oo, solve 


min{ 9 r 2 : | x(e r , 0, 0, #)} = q T x r+ . 




min {q T x | x[(£™ m ~ lr)e r , 0, a r ,/3 r ]} = q T x r+ (-^ ,n + | r )- 
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Else (let = oo if f3 r (i) = -boo, /?*(i) = 0 otherwise) and solve 



mm{q T x\ x{-er,0,0,p[)} = -q T x r . 

If Step 1 was entered with x x± = (0, ... ,0) for all i, go to Step 2; otherwise, go to Step 4. 



Step 2: For t = 1, . . . , n, solve 




« r+1 (0 


= max — x°(t) — ar ,+ (*)y J+ - ^"(O^ 

&r+ 1 JT^r+l 


subject to 


y y+ -y y - =£;-£; 



P r+l (i) = min- z°(t) - ^2 x J+ (i)y 3+ - ** (Oy 7 + 00 



subject to 


jVr+l jVr+l 

y y+ - y y ~ =£/-!; 

ef n < e* *r+i. 



Step 8: If r < m, let r = r + 1 and go to Step 1. 
Otherwise, go to Step 4. 

Step 4-‘ Find 



a*(i 


mi mi 

0 = max - z°(0 - JZ * /+ (*V + “ ^"(Ov* - 

3=1 3=1 


subject to 


y ,+ - y y ~ =6-£y 

£f n <^<^ ax -y = i,...,m 1) 


/no = 


mi mi 

= min - z ° (0 - XXn0y J+ +c ( t ) 

3=1 3 = 1 


subject to 


y y+ - y y ~ = €;-£; 

£f n < 6- <c>n n--> m i> 



for i = 1, . . . , n. 



Let x* = argmin {q T x | x(0, <£ max , a*, ft*)}- Find a conformal realization of x* (Rockafellar [1984, p. 455]), so 
that 

x * = a k x k a k > 0 , 

such that x*{i) > 0 => x k (i) > 0 and x*(z) < 0 => x k (i) < 0, and x*(z) = 0 => x k (i) = 0. An algorithm for 
finding such a realization is the “painted index algorithm” in Rockafellar [1984, p.476]. Paint all columns Aj 
of A such that 

{ white if x*(j) > 0, 
black if x*(j) < 0, 
red if x*(j) = 0. 

Let k — 1. Pivot until a Tucker- tableau is reached in which there is a compatible column. This will always 
be possible in our case. Let the compatible column be Ay, and let F be the set of indices for the basic 
columns in the final Tucker tableau. We now have that 



£>•4(0 + Ay = o. 

ieF 



If Aj is white, let 



(A'.(i) if ieF, 
z it( l ) = S 1 if * = j, 

1 0 otherwise. 

If Aj is black, reverse all signs in x* k . (Note that the sign convention in a Tucker tableau is opposite of the 
convention in the standard simplex tableau.) 

Let ak = min{x*(i)/x£(i),x£(z) ^ 0}, x*(i) = x*(i) — akX k (i) and re-paint every column for which x*(i) = 0 
red. 

If x* 0, let k = k + 1 and repeat. Otherwise, go to Step 5 with the conformal realization a k x l- 



Step 5: Using the cost coefficients q T x i± i find This amounts to performing mi simple line 

integrals. 



Step 6: If x* (i) > 0 (so that x k (i) > 0, Vfc), we are using a variable x(t) with random capacity /3* [i) +<}>{(> 

If x*(i) < 0, we are using a variable x(z) with deterministic capacity a*(t)(< 0). We shall in the following 
assume that each variable x(z),such that x*(z) ^ 0, has associated with it a random arc capacity <f>*. 
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If x* (i) < 0, we have Pr{<^ = a*(i)} = 1, if x*(z) > 0, <p* = <f> + /?*(i). For each k = 1 let 

q k = £g(i)x£(i)(< 0). Sort the primal supports x k such that qi < q 2 < • • ■ < <7jt. Let k = 1, /? = 0 (where 
p will become Eff(<£)). 

Step 7: Let P = {i | x£(i) ^ 0}. Consider the random variable 

Pk = max{0,min P {^/xfc(i)}},^ fc € [0,a fc ]. 

Find (This work amounts to increasing the capacity of each conformal flow until the first variable 

capacity is met. This continues on each conformal flow. Details are given for the network case in Wallace 
[1987b].) Let p = p + E qkfik, and <f>* = <j>* — o^x £. If A; = K or if q^+i = 0, stop with p = E otherwise 
let A; = A; + 1 and repeat Step 7. 

End. 

The value obtained in Algorithm 1 is indeed an upper bound on the expected linear program value. 
Theorem. The value SPLU= E <f>)} obtained in Algorithm 1 is an upper bound on Q = E^[<3(£, <£)]. 

Proof: The proof requires only showing that x = x° + ^(x J + (£y — £) + + x J ~(£ — £y) + ) + £) is 

feasible in x(£> <£*> c )- This is obtained by noting that the definitions of x r± , a r , and fi r in Steps 0 to 2 
maintain feasibility for <f>* .1 

The algorithm as described above is our basic version. We prove certain properties of it in the next 
section. In Section 5, we present alternative versions of some of the steps in Algorithm 1. 

3. Properties of the Upper Bound 

The purpose of this section is to show that the upper bound presented in this paper has some desirable 
properties and to relate the procedure to other bounding methods. 

3.1 Exact Bounds for Linear Problems 

All other bounds used in stochastic programming are exact whenever Q(£, <£) is linear in £ and <j> over 
the support of the random variables. This is true of the Madansky upper bound, the piecewise linear upper 
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bound in the pure network case (Wallace [1987b]), the linear upper bound on the expected max flow in a 
network (Wallace [1987a]), the Jensen lower bound, and the sublinear approximation in Birge and Wets (if 
the random variables have unbounded support). All acceptable bounds should have this property. 

Property 1: The bound SPLU given by Algorithm 1 is exact if Q(f, <f>) is a linear function. 

Proof: Assume Q{l- } <f>) is linear in £ and </> and that the reduced cost of a non-basic variable is always 
different from zero (a dual non- degeneracy assumption). Then EJ7(£,<£) = EQ(£,0). Of course, if Q is 
linear, it can be written as 

mi 

Q{t, 4>) = Q{1 o) + E fkitk ~ &) + E M/- 

k=l j 

Clearly, Step 0 provides us with <2(£, 0). Also, if Q is linear, a 1 4- e < 0, fi l -f / > 0 in Step 0, since the basis 
corresponding to Q(£, 0) is feasible for all £ £ 5. Hence, / % = qx = —qx k ~. Therefore, if Q is linear, the 
algorithm will discover the coefficients of £ in Step 0 and then go to Step 4. 

Let us define a variable i to be stochastic if <£™ ax > 0, otherwise, it is deterministic. Consider the 
conformal realization of x* = ^ockX * h . First note that x* k is an elementary vector (Rockafellar [1984, p.453]). 
This means that there is no way to split x k into two or more other vectors where at least one has fewer 
non-zeroes than x k . 

Assume there exists an elementary vector y such that y(i) ^ 0 for more than one stochastic random 
variable. Then fix the value of fa at 0 for all variables except for those with y(i) ^ 0. Then, Q would not be 
linear. (Compare with the random variable fi in Step 7.) Hence, if Q is linear, there is no elementary vector 
with more than one stochastic variable. 

Now, assume that we have found two elementary vectors t/i and t /2 , such that they share the stochastic 
variable % (i.e., t/i(t) ^ 0,y 2 (i) ^ 0). Also assume that q\/y\(i) ^ < 72 / 3 / 2 (0* (The variable qi defined as in 
Step 6 .) Let all <j>j — 0 for i 7 ^ j. Then Q is not linear in variable t, because the marginal gain of increasing 
<f>i is not the same in both elementary vectors. Hence, two elementary vectors can only share a stochastic 
variable if qi/yi(i) = < 72 / 3 / 2 ( 0 * (This corresponds to two circuits in a pure network that have the same cost 
and share an arc with a random capacity.) Of course, hi = <7i/yi(0* 

Hence, if Q is linear, no elementary vector x* k has more than one stochastic variable and two elementary 
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vectors can only share a stochastic variable if they have the same cost (in the sense described above). Since 
Step 6 only creates element ary vectors x£, the random variable fik in Step 7 is linear in its single random 
variable. Hence, our method produces the exact solution.^ 

3.2 The Bound is Polynomial 

The Edmundson-Madansky bound requires that Q(£, <f>) be solved in all extreme cases of £ and <f. There 
are 2 mi "*' ni such points; hence, the method is exponential in the number of stochastic variables. Only for 
very moderate values of n x and m x is it possible to apply this bound. 

The major goal of this paper is therefore to find a good upper bound that can be computed in a number 
of operations that is polynomial rather than exponential in the number of random variables. 

Property 2: Algorithm 1 calculates SPLU in a number of operations that is polynomial in the number of 
random variables. 

Proof: The amount of work is in the worst case: 

Step 0: 1 LP (a 1 ,/? 1 can be found by inspection). 

Step 1: 2m x LPs. 

Step 4: 1 LP to find x* . The conformal realization is independent of n x and mi. (The worst case is n LP’s, 
n < n x .) 

Step 5: The integration is a constant amount of work for each random variable. 

Step 7: Finding amounts to checking the mi * max{m*} (in the worst case) possible values of The 

value mi is the total number of possible values for fc. This has to be done not more than n times (since the 
number of zeroes increases by one for each A;). 

Hence, the algorithm is linear in n x and miJ 

3.3 Relation to Networks 

The method presented in this paper is closely related to the network method in Wallace [1987b]. The 
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major difference is in Step 1, where we only solve two networks in the network case and not 2 as here. 

Below is a short network interpretation of some of the vectors and scalars used in the algorithm to help 
in its understanding. 

Step 0: The variable rc t+ shows how the flow changes on the basic arcs as the supply at node i is increased 

by one unit or (the demand is decreased). Hence, 

( +1 if arc j is a forward arc on the path from node i to the slack node, 

— 1 if arc j is a reverse arc, 

0 if arc j is not on the path. 

x x ~ is similarly defined for increased demand (or decreased supply). 

a 1 ^’) > 0 implies that with the chosen set of paths (rr^) there are supply/demand combinations that give a 
negative flow on arc i , even when we disregard node 1. 

/^(i) > 0 implies that with the chosen paths, there are supply/demand combinations that overuse arc i even 
when we do not consider node 1. 

Step 1: xf are still paths, but not along a basis. Both basic and non-basic arcs are used. If Step 1 finishes 
successfully, we have actually replaced the original network by a star-shaped network (where the slack node 
is in the center of the star). The arc going from the center node to node i has unit cost the arc in 

the other direction has unit cost qxf . The way we have used a r and has guaranteed that whatever 
combination we get of supply and demand, sending that flow along the paths xf would be feasible and cost 
the same as in the star-shaped network. 

Hence, we have found an upper bounding simple recourse problem (Wets [1983]) . In stochastic programming 
this approximation depends on the actual value of the first stage decisions (as in the recourse function in 
the introduction). Hence, in some sense, it is a local approximation. 

Step 4: a* < 0 shows how much flow can be sent along the original arcs in the negative direction without 
making that total flow negative (whatever the supply/demand is). Similarly, ft* shows how much is left of 
the capacity in the arcs in the worst case. 

x* is just a circulation in the network, and x* k are circuits of minimal length (in terms of the number of 
arcs in them), shows how much flow the circuit can take (or, more precisely, how much flow it has been 
allotted.) 
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Step 7: Pk is a random variable describing the capacity of circuit k< 



3.4 Relation to Sublinear Approximations 

The separable piecewise linear upper bound is also a generalization of the ray function approximation in 
Birge and Wets [1986a] and its extension in the sublinear approximation in Birge and Wets [1986b]. These 
procedures find the value of Q(£, <f>) in different coordinate directions to again obtain a separable function 
that can easily be integrated. The approach in Birge and Wets [1986b] uses varying choices of the coordinate 
system that leads to an extension of the SPLU bound given here. This extension would involve solving for 
x J + and x°~ in different directions so that a variety of bounds could be obtained. 

The ray function approximation amounts to solving for 

q T x 3 * = min {q T x | Ax = e° , x > 0} 

and 

q T x J ~ = min{g T x | Ax = —t\x > 0}. 

These values of x J+ and x J ~ are then used in J7(£, <f>) as in SPLU. The extension is to use the elements of 
other coordinate systems in place of in the definitions (i.e., use some vectors d? that form a basis for 
$l n ). This procedure can be used in Algorithm 1 to obtain an alternative bound. 

The sublinear approximation with varying directions has been found to produce accurate approximations 
in a variety of examples. The advantage of the SPLU bound is that it applies to bounded regions so it may 
be used on partitions of the support of the random variable in a refinement procedure in solving a stochastic 
program. Algorithm 1 also incorporates the procedures for handling random bounds that often arise in 
practical examples. 

3.5 Finiteness 

There is no guarantee that our upper bound is finite, i.e., that all linear programs that must be solved 
are feasible. An infinite bound of course results if EQ(£, <f>} = -boo, i.e. the problem itself is infeasible, but 
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it can also be that E<3(?,^) < +oo, whereas Ef/(£,$) = -foo. This is not always avoidable. We note that 
the only other polynomial upper bound, the ray approximation, is never better than our bound (assuming 
the possible extensions mentioned above), and that exponential bounds may be necessary in some cases. 

3.6 Partitioning 

When approximations, such as the one in this paper, are used in two-stage stochastic programming, a 
comparison is made with a lower bound (EL) usually based on Jensen’s inequality. Then, if E U - EL is 
too large (according to some rule), the support rectangle (for independent random variables) is partitioned 
into smaller rectangles called ce//s, and the bounding procedures are applied to these cells, which in turn are 
weighted by their probability. 

Hence, whenever a partition is called for, one must decide which cell to partition and, along which 
coordinate direction, to perform the partition. With an upper bounding method that can take on the value 
+oo even for a feasible problem, one should clearly partition the cell where E U = +oo, along the coordinate 
direction that was being treated when the infeasibility was discovered. This provides a dynamic scheme in 
which the algorithm is applied on each cell until either an infinite value is obtained or the difference between 
lower and upper bounds is above the acceptable threshold. A partition is made in either instance. Partition 
strategies are discussed in Birge and Wets [1986a], Birge and Wallace [1986] and Frauendorfer and Kali 
[1986]. 

4. Examples 

In this section, we first present a small example to illustrate the bound. We then give computational 
results on a larger problem from energy modeling (Louveaux[l987]). 

4.1 A Problem with Two Random Variables 

The first example is a problem with two random variables and without random capacities. We wish to 
find bounds on E Q(£) where 

Q(£) = min x\ + x 2 £3 + z 4 + IOX 5 + 10x6 (4.1) 



13 



subject to 

xi +3x 2 +*3 -x 6 = £1 (4.2) 

3xi +X2 +X4 — X6 = £2 (4.3) 

Xx,...,x 6 > 0, 

where £1 and £2 are uniformly distributed on [1,4]. This problem is illustrated in Figure 1, where A, refers 
to the ith column in the constraint matrix of (4. 2-4. 3). We follow Algorithm 1 step by step. 

Step 0: 

(i) Find Q(£) = 1.25. 

x° = (0.625,0.625,0,0,0,0). 

x 1+ = (—0.125,0.375,0,0,0, 0); x 1 ” = (0.125,-0.375,0,0,0,0). 
x 2 + = (0.375, -0.125, 0,0,0,0);x 2 “ = (-0.375,0.125,0,0,0,0). 
a x (l) = -0.0625; a 1 (2) = -0.4375;(Note : $ r [i) = + 00 .) 

Now e(l) = 0.1875 > -a 1 (l) = 0.0625, so go to Step 1. Note in Figure 1 that we have essentially moved 
along the vertical line through £. The bound a 1 recorded the (negative) minimum multiples of the vectors 
A. 1 and A .2 for points along that line. The value e(l) recorded the greatest change in the multiple of A.i 
from the multiple for £ for other points along the horizontal line through £. The function is not linear 
because this change is greater than the minimal multiple (ai(l)) for movement in the vertical direction. 

Step 1 . 

Solve min {q T x | x[1.5«i, 0, (-0.0625, -0.4375, 0, 0,0,0), 00 )]} = 1.125 = (0.75) * (1.5) = g T x 1+ (^ ax - £ 1 ), 
where x 1+ = (-0.0625,0.1875, 1.0, 0,0,0), and 

min{g T x | *[-1.5ei, 0, (-0.0625, -0.4375,0,0, 0, 0), 00 )]} = 1.375 = (0.9167) * (1.5) = g r x 1 ”(6 - £5 nin ), 
where x 1- = (-0.0625,-0.4375,0,0.625,0.125,0). Next go to Step 4. 

Step 4: 

We can skip this since there are no random bounds. Step 5 then is the terminal step. 

Step 5: 
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Here we compute Ef/(f) = 



<?(£)+/ * I+ (6-£i)<tf , (&)+ / - liWli) 

J £i>(i J£i<£i 

+ [ Z 2+ (f2 - |2)<iF(^2) + [ Z 2_ (|2 — 

J •'£a<^3 

= 1.25 + 0.5(0.75)(0.75) + (0.5)(0.9167) (0.75) + 0.5(0.25) (0.75) + (0.5)(-0.25)(0.75) = 1.875. 

End. 

So, we have SPLU = 1.875. We compare this with the Edmundson-Madansky (EM) bound. In this 
example, the EM bound assigns equal weights to the values of <?(£) at each of the extreme points of 5. 
Hence, 

EM = 0.25 * (<?( 1 , 1) + Q( 1 , 4) + Q{ 4, 1) + <?(4, 4)) = 1.625. 

The EM bound is better than the SPLU bound but this difference may be eliminated by refinements of the 
SPLU bound. We describe possible refinements in Section 5. 

4.2 Computational Results for an Energy Model 

The usefulness of the SPLU bound is best demonstrated on a practical example in which the number of 
random variables varies. We wish specifically to observe the performance of SPLU relative to the EM bound 
as the number of random variables increases. The performance is measured in the sharpness of the bound 
and the computational effort. As a practical example, we consider the small energy model in Louveaux 
[1987]. We do not consider random bounds because that is directly analogous to the network case discussed 
in Wallace [1987b]. 

In this example, we have four technologies which can be used to satisfy three demands at varying costs. 
High cost “backstop” technologies are also available to satisfy demand so the problem is feasible for any 
demand realization. The randomness occurs in the capacity of the technologies and the demands. This allows 
from one to seven random variables. The examples were also chosen with varying ranges (narrow, medium, 
and wide) on the random variables resulting in twenty-one sets of examples. We assume uniform distributions. 
This assumption favors bounds (such as the Edmundson-Madansky bound) that place weights at extreme 
values since other distributions generally have more mass around the center of the support. 
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The experiments were conducted on the Amdahl 5860 at The University of Michigan Computing Center. 
The SPLU and EM bounds were both implemented in FORTRAN codes using the same linear programming 
routine LPM-1 (Pfefferkorn and Tomlin [1976]). Each bound was computed for each of the twenty-one test 
problems. The Jensen inequality lower bound was also computed to determine the values of the upper 
bounds relative to the lower bounds. The results are given in Table 1. 

The results in Table 1 show that the polynomial bound SPLU does not generally provide as accurate a 
bound as the EM bound, but that as the number of random variables increases the computational time in 
SPLU increases much less rapidly than the time for EM. In these examples, the growth of time for SPLU 
is indeed approximately linear (gaining ten milliseconds for each random variable), while the time for EM 
approximately doubles as each new random variable is introduced. This demonstrates that the real advantage 
of the SPLU bound is in problems with a large number of random variables where the EM bound cannot be 
computed. These results are comparable with the results in Wallace [ 1987b] for networks. 

Refinements are also possible to reduce the error in SPLU. In Section 5, a refinement scheme using 
parametric linear programming is introduced. The use of different coordinate directions is another possibility 
as mentioned above. For unbounded ranges, the resulting sublinear approximation values were reduced up 
to thirty percent from the coordinate direction values for a similar set of test problems (Birge and Wets 
[1986b]). 
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PROBLEM 1 


TIME 2 




VALUE 






EM 


SPLU 


EM 


SPLU 


Jensen 


l-NAR 


9 


10 


182.75 


182.75 


182.75 


1-MED 


10 


18 


220.50 


220.25 


220.00 


1-WID 


10 


21 


385.50 


341.50 


297.50 


2-NAR 


16 


20 


183.38 


183.06 


182.75 


2-MED 


17 


26 


220.50 


220.50 


220.00 


2-WID 


22 


35 


389.85 


389.10 


297.50 


3-NAR 


22 


28 


183.38 


183.38 


182.75 


3-MED 


25 


38 


221.38 


222.50 


220.00 


3-WID 


33 


40 


433.60 


439.58 


297.50 


4-NAR 


51 


41 


184.09 


185.50 


182.75 


4- MED 


44 


41 


227.22 


255.18 


220.00 


4-WID 


47 


45 


434.26 


469.50 


297.50 


5-NAR 


94 


52 


184.19 


186.44 


182.75 


5-MED 


87 


51 


227.41 


278.38 


220.00 


5-WID 


75 


54 


434.35 


499.18 


297.50 


6-NAR 


163 


54 


185.58 


192.49 


182.75 


6- MED 


193 


61 


235.91 


303.35 


220.00 


6-WID 


149 


72 


443.52 


524.30 


297.50 


7-NAR 


329 


64 


186.23 


215.48 


182.75 


7-MED 


366 


70 


236.27 


328.35 


220.00 


7-WID 


304 


76 


444.12 


549.30 


297.50 



1 - Number of random variables - range of random variables. 

2 - CPU milliseconds. 

Table 1. Results for Edmundson-Madansky and SPLU bounds. 



17 



5. Extensions and Conclusions 



The SPLU bound can be refined in a variety ways. The use of other coordinate directions may be 
possible, but it is best used when linear transformations of the random variables have a known distributional 
form as is the case for normally distributed random variables. As mentioned above, a common procedure 
is to partition the support of the random variables and to apply the bound on each of the partitions. Here 
we give a parametric programming approach that can obtain more accurate results without partitioning the 
random variables. The following modifications of Algorithm 1 provide this basic bound. 

Algorithm 2 

Substitute the following steps into Algorithm 1 to obtain 

Vti,4>)-Q(lo) + H(4)+ E -(+)$)■ 

€.>(<)! 

Step 1 ' .Solve the parametric linear program 

min{< 7 r a; | x(« r) 0, a r ,/9 r )} 

for c G [0, — £ r ] (or € G [0, oo) if £ r is unbounded). This generates a piecewise linear function /+(ee r ) 

with break points {0, €i, . . . , €r}, and with slope values, q T x J"* - , . . . , q T x^ + . 

Then, solve the parametric linear program 

min{q T x | 0, a r , £ r )} 

for e G [0, £ r — £™ lu ] (or c G [0, oo) if £ r unbounded.) We then obtain a piecewise linear function /“ with 
breakpoints {0, c 1 , . . . , c T }, and with slope values, q T x r {~ , . . . , q T XjT . 

Step g and Step : Substitute mint{xt + (z)} for x J+ (z) and mint{x/~ (z)} for x 3 ~(i) in the definitions of 
a r+1 (z) and a*(z) and substitute maxt{x^ + (i)} for x J+ (z) and max t {x^”(z)} for x J ~(z) in the definitions of 
P r+i {i) and 0*(t). 

The changes in Step 1' of Algorithm 2 lead to better bounds if a r and f3 r are the same and q T x[ < q T x r T . 
In Algorithm 2, the approximation obtains as low a value as possible for all for changes in the rth direction 
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given the values found for movement in previous directions. In Algorithm 1, the approximation just uses the 
extreme values of £ r . 

This difference can be seen in the example from Section 4.1 which is illustrated in Figure 2. The dashed 
line corresponds to the function used in Algorithm 1 by using the extreme values. The solid line corresponds 
to the functions /+ and /“. The new bound is 

SPLU' = E[{7'(£, <£)] = 1.449. 

We note that SPLU' is now below the EM bound value of 1.625. 

The bound from Algorithm 2 is not always better than SPLU because the bounds may change for 
different values of r, i.e. a r+1 may increase and /? r+1 may decrease. Although this difference appears to rarely 
make SPLU' worse than SPLU according to our limited conmputational experience, it may be advantageous 
to guarantee that a bound at least as good as SPLU is obtained. This guarantee is accomplished in the 
following modification of Step 1'. 

Step 1 " . Solve 

mm{g r x|x(U r max -|rK,0,a r ^ r } 



Let 



a J(i) r *(*') if x < 0 , 

1 0 otherwise. 



and 



B r U) = [ *(*) 1 > 0 

l 0 otherwise. 



Then solve the parametric linear program 



min{g T x | x((«v)e r ,0, a',/?!)} 



to obtain f?(te r ) as explained in Step 1'. 



19 



This modification of Algorithm 2 results in bounds that are at least as sharp as SPLU and can still 
benefit from the parametric program as in the example given above. The key benefit of the SPLU bound 
that the computational effort only grows polynomially with increases in the number of random variables is 
maintained. We have demonstrated how this improvement results in reduced times on one set of examples 
and that the greatest value of the SPLU bound may be in cases where the EM and other exponential bounds 
cannot be reasonably computed. The refinements mentioned above may allow the SPLU bound to be even 
more useful in the solution of practical stochastic linear programming problems. 
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Figure 1 . Example with two random variables. 
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Figure 2. Parametric linear program bound. 
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